Self-organized criticality in the Kardar-Parisi- Zhang-equation 
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Kardar-Parisi-Zhang interface depinning with quenched noise is studied in an ensemble that 
leads to self-organized criticality in the quenched Edwards-Wilkinson (QEW) universality class and 
related sandpile models. An interface is pinned at the boundaries, and a slowly increasing external 
drive is added to compensate for the pinning. The ensuing interface behavior describes the integrated 
toppling activity history of a QKPZ cellular automaton. The avalanche picture consists of several 
phases depending on the relative importance of the terms in the interface equation. The SOC state is 
more complicated than in the QEW case and it is not related to the properties of the bulk depinning 
transition. 
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One of the more interesting developments in non- 
equilibrium statistical physics has been the realization 
that inhomogeneous systems give rise to new, rich behav- 
ior compared to those in which translational invariance is 
present. A prime example thereof is the criticality found 
in so-called sandpile models, coined 'self-organized' due 
to the absence of an obvious control parameter . 

Usual 'self-organized critical' SOC sandpiles have two 
central features. First, the system is driven very slowly 
so that individual avalanches occur on a separate, fast 
time-scale. Second, the drive (like the addition of grains) 
and the dissipation (loss of grains) occur non-uniformly. 
Typically the system is forced, on the average, in an uni- 
form fashion but open boundaries create dissipation, in 
sandpiles via the loss of grains. These two characteristics 
are the keys to describe, and the earmarks of SOC (for 
early work see ||) [||^]. 

Essential to a wide class of SOC models is that the 
movement of grains has a Laplacian character. Con- 
sider now the history of activity, topplings in such mod- 
els. The activation of a site means that it has re- 
ceived enough grains, z(x,t), from outside or through 
redistribution, to overcome a local activation threshold: 
z(x,t) > z c (x,t). For a large class of models, this im- 
plies that if H(x, t) defines the time-integrated activity, 
H = f p(x, t)cLt = f Q(z(x, t) — z c (x, r))dr, then after 
discretization H obeys the discrete interface equation 



H(x,t + 1) 



H(x,t) + 1, f(x,t) > 0, 
H(x,t), f(x,t)<0, 



(1) 



where / denotes the local 'force' acting on the interface 
W. This can be recast in the form 



AH 

~aT 



(vV 2 H + r]{x, H) + F{x, t) - F c 



(2) 



The Laplacian comes from the average redistribution of 
grains, and a constant v takes into account the specific 
updating rules of the system. r\ is a quenched noise 
term that originates from the random rules of sandpiles. 



F(x, t) measures the external drive and changes on the 
slow-time scale, and F c is the toppling threshold writ- 
ten as a force acting on the interface. By casting the 
dynamics of a sandpile in the form of a driven inter- 
face in a random enviroment (rj) one sees that the SOC 
critical point is reached under the combined action of a 
slowly-increasing external force F(x, t) and the boundary 
conditions, H = 0, which amount to the loss of grains. 
This is a slightly unusual case of a depinning transition 
of the quenched Edwards- Wilkinson universality class, 
since the macroscopic shape of the interface is parabolic 
(in ID) The avalanche picture in an ordinary depin- 
ning transition J| |8| is however different since the trans- 
lational invariance is here broken as manifested already 
in the average interface profile - which tells us that the 
activity density peaks in the bulk, in the center of the 
system. 

The connection of SOC to non-equilibrium phase tran- 
sitions can be made both in the interface context and 
in the realm of absorbing state phase transitions, ex- 
emplified by the contact process and various reaction- 
diffusion models One way to classify SOC models 
is to consider the noise correlator in the presentation at 
hand. Such analysis however assumes certain symme- 
tries included in the representation (the QEW Langevin 
equation; the Langevin equation describing the active 
site density). A pertinent question thus arises: what 
happens if the ensemble, in which SOC is observed, is 
varied with respect to its fundamental symmetries? In 
this work we illustrate this by analyzing the outcome for 
the Kardar-Parisi-Zhang universality class. The differ- 
ence between the QKPZ and the EW equations (or a 
non-linear Langevin equation compared to a linear diffu- 
sion equation) is that the diffusion equation is insensitive 
to most initial profiles due to the fundamental symme- 
tries of the equation, while in the QKPZ case the growth 
component perpendicular to the interface slope is always 
non-zero. This means that there will be a new history 
dependence in the QKPZ interface which is absent from 
the EW dynamics and from usual absorbing state phase 
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transitions. 

Consider, in analogy to the sandpile interface repre- 
sentation, the quenched KPZ equation which reads 

A JT \ 

— = Q{W 2 H + -{VHf + „(x, H) + F(x, t)), (3) 

where the nonlinearity proportional to A has been added 
to the QEW equation. Note the step function that en- 
forces the constraint that the interface cannot propagate 
backwards. In the case of the sandpile equation, it also 
discretizes the interface velocity to zero or unity which in 
the absence of mixing by randomness creates an impor- 
tant Abelian symmetry and thus becomes crucial in that 
case jDJ. The external drive is in the following taken 
to be uniform in space, F(x, t) — > Fit) and the noise to 
be quenched and delta-correlated, (rj(x, H)i~i(x! , H')) = 
2DS d (x - x')5(H - H'). Equation @ has to be aug- 
mented with suitable boundary conditions in order to 
reach a steady-state. As in the QEW or sandpile ensem- 
ble, we use the selection H(x, t) = at x = 0, x = L + 1 
(in 1+1 dimensions) instead of e.g. periodic boundary 
conditions. In the case of such ordinary boundary con- 
ditions Eq. (H) has a depinning transition at a critical 
F c . The exponents depend on whether A approaches zero 
with AF = F — F c — » or not, and in the latter case we 
have a true QPKZ depinning transition fl(i| , |l2]| . 

An interface evolving from a flat boundary condition 
will develop a non-zero roughness and finally stops due 
to the influence of boundary conditions. The roughening, 
the development of average local slopes, depends on the 
sign of A. To reach SOC-like conditions, the separation of 
timescales, or a vanishingly small but positive interface 
velocity (average activity level) v = + one has to choose 
the right way of driving the system, F(t). For F(t) > 
the interface becomes unstable if A > and ^ > 0. This 
breaking of up/down symmetry, absent in the continuum 
QEW model, shows the first sign that the possible SOC 
state in the QKPZ equation is different. If A < it be- 
comes posssible to take the SOC limit of vanishingly slow 
external drive, 

F(jt)=Ct, C^O (4) 

from which follows that a separation of time-scales hap- 
pens, similarly to the QEW or normal sandpile case. In 
the following we investigate the ensuing interface dynam- 
ics. 

To this end, a cellular automaton jl3) version of Eq. (||) 
is defined as follows. Let fi be a microscopic force at the 
discretized sites x —> i, such that for fi > i topples, fi 
reads in ID 

Mt) = v [H i+1 (t) + Hi-x(t) - 2H.it)] 

+ ^ ([fli(t) - Hi-i it)} 2 + [Hi+t® - Hiit)} 2 ) 

+Fit)+r HtHi . (5) 

This can be interpreted, in analogy to SOC ones, as 
a QKPZ sandpile where a site gives 2v particles or 



units of energy to its neighbors once it topples, and 
the toppling criterion also compares the integrated ac- 
tivity at site i to its neighbors. In the interface pic- 
ture Hi(t + 1) = Hi(t) + Q(fi(t)). With this definition, 
Hix,t) — f pix,r)dr maps the discrete QKPZ equa- 
tion to the 'sandpile' and vice versa. Notice that the 
nonlinearity, proportional to A, above is discretized in a 
particular way, to avoid the numerical instabilities of the 
QKPZ equation in a reasonable fashion p4). 
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FIG. 1. Typical interface configurations for a ID system of 
500 sites. The one with the flat top shows an unstable one 
from the initial growth event before the first 'triangle' is cre- 
ated. The other configurations are stable, from every 20000th 
time step. 



Fig. shows a series of snapshots from a simulation 
of the QKPZ equation with r)i Hi — i<? with equal prob- 
ability. This describes random toppling thresholds that 
change after each toppling in the sandpile language. The 
system is driven in such a way that F(t) is increased af- 
ter an avalanche has relaxed (defined in the usual way 
by the stoppage of all activity) , so that one site becomes 
active, and the interface advances. It is evident that the 
dynamics of such a system consists of at least of the devel- 
opment of an initial profile, and the subsequent phase(s). 
The macroscopic shape of the QKPZ interface is to first 
order triangular in 1+1-D as expected with A < 0. The 
slope follows from the different terms in Eq. (|J). For a 
flat slope the surface tension is negligible, and the slope 
arises from the balance of the nonlinearity and the driv- 
ing force. This gives (| V/i |) = \/2(t/X. The development 
of such a profile, determined by the symmetries of the in- 
terface equation, is analogous to the development of the 
steady-state grain configuration in ordinary SOC sand- 
piles from an empty system. Next we discuss the phase 
that appears after the initial transient. 

As in ordinary SOC, avalanches are characterized by 
the duration T, area I, and number of topplings s (the 
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number of times sites advance). The initial growth, giv- 
ing rise to the triangular shape, is neglected. We obtain 
for the avalanche sizes in this regime the result shown 
in Fig. |^ for four system sizes. The size distribution has 
asymptotic power law decay with a measured exponent 
of r s = 2.5 ± 0.1 with P(s) ~ s~ Ts . Due to the large t s 
exponent the average avalanche size is not dependent on 
L, unlike in usual SOC sandpile models. 
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FIG. 2. Probability distribution function of the avalanche 
sizes for ID ensembles. System sizes are 250, 500, 750, and 
1000 for A, o, □, and O, respectively. The parameters are 
v — 20, A = —4, g — 10. Dynamics was only followed for 
configurations with average interface slopes less than 10 (see 
text). The parallel line indicates the region of power law fit- 
ting resulting in an exponent of —2.5 ± 0.1. 



However, the avalanches are time-translation invariant 
only in the average sense, due to the role of the average 
interface slope to. The avalanche sizes are plotted ver- 
sus time on Fig. |. On the same diagram, the standard 
deviation D of the local slope is also shown. 
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FIG. 3. Avalanche sizes s versus time and the standard 
deviation of the local slope D. v = 20, A = —4, g = 10 
and the system size is 1000. The correlation between the two 
quantities is obvious. 



The largest avalanches take place when the roughness 
of the interface D is the largest. Smaller, but still rel- 
atively large avalanches develop where the sides of the 
triangular interface are smooth. The local slopes, i. e. 
\Hi — Hi_i \ alternate between n and n + 1, where n is a 
positive integer. It is not surprising that a smooth inter- 
face evokes large avalanches, because there are no fluc- 
tuations along the surface that would stop an upwards 
propagation of the avalanche. When D is the small- 
est, the first few avalanches start from either bottom of 
the triangle and propagate up to the top, increasing the 
interface height by one at every site. The really large 
avalanches occur at the maxima of D. They typically do 
not increase the average gradient of the side, but only the 
height of interface evenly along the side. In spite of the 
extra slope-dependence, there is no characteristic size as 
is typical of SOC and the correlation length equals the 
system size. The period of the avalanches is simply re- 
lated to the squared system size. This is since the area 
covered by all the avalanches between two integer values 
of m is proportional to L 2 , and, the average avalanche 
size and duration are of the order of unity: (s) ~ 0(1) 
and (T) - 0(1). 

The support or the area of the avalanche is as usual 
defined by the number of sites that have toppled at least 
once during an avalanche. Notice that an avalanche al- 
ways tends to spread upwards, and thus in ID the start- 
ing and final, extremal sites usually give the area directly. 
Sometimes, the Laplacian makes it so that an avalanche 
can propagate over the top of the 'triangle'. We find 
that the probability of starting an avalanche linearly in- 
creases with the location of a site and the probability 
of an avalanche to stop grows linearly also as the top is 
approached. 

The area of the avalanches, defined as I = i r — %i + 1, is 
expected to follow a similar power law as the avalanche 
sizes do and one obtains that 77 ~ 3.0 ±0.1. The 
avalanche lifetime T is the time that elapses between two 
successive stationary states in the system and also be- 
haves as a power-law, P(T) ~ t~ 29±0 - 2 . Due to the fact 
that the average avalanche properties are defined by the 
small avalanches no analogy of the usual scaling relations 
between the exponents seems to exist. 

The above results have, however, been obtained in 
just a window of time. Consider an interface (side of 
the triangle) with an uniform slope m and an avalanche 
which has just started to propagate and reached site i. 
The local force at site i, at time t, becomes fi(t) = 
v(m - 1 - to) + I [(to - l) 2 + m 2 ] + F + n t 2 = -v + 

^(2m - 1) - ^2to 2 + F + % and at i + 1 f i+ i(t) = 
v[m - (m - 1)] + I [to 2 + (to - l) 2 ] + F + n l+1 = 
v + ^(2m — 1) — Jj-2m 2 + F + r/i+i- Assuming that on 
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average on an even interface F balances the local forces 
for slope m (F w ^-2m 2 ) and neglecting the noise fluc- 
tuations, we find a critical slope 



2// 



(6) 



This means that the cell at i may stay active for 
the next updating step at time t + 1 when m > m cr . 
The natural progress of avalanches is now such that the 
site at i + 1 is normally always activated if the inter- 
face is smooth enough. In this phase the site at i + 1 
gets active and topples at t + 1, but at the same time 
Hi — > i/i + 1 as well. Thus, a short region of the inter- 
face with a smaller gradient will be formed and starts to 
propagate upwards, giving rise to very large avalanches. 
Since these are more difficult to induce, the majority of 
the avalanches are of size unity, and tend to smoothcn 
out the interface fluctuations. The avalanche statistics 
data presented previously is restricted to systems with 
a maximum slope of 10: the critical slope for those sys- 
tems is m cr « 2 • 20/4+ 1/2 = 10.5. In this regime above 
the critical slope m cr , the scaling structure of avalanche 
behavior is lost, and the growth process produces either 
extremely large or very small avalanches. This feature 
makes it prohibitively difficult to take this regime un- 
der closer scrutiny with numerical simulations. The time 
to reach m cr can be estimated in two ways depending 
on whether one uses macroscopic time or 'SOC time 
(counting the number of avalanches). With a uniformly, 
but very slowly increasing force F the critical slope is 
reached at the same time regardless of the system size 
L. In contrast, since the average size of an avalanche is 
independent of L (recall that t s > 2), it follows that the 
number of avalanches n cr to reach m cr will depend on L 

aS Tl(yr L . 

The previous results illustrate the nature of SOC in 
the QKPZ universality class. The combination of slow 
drive and boundary pinning chooses a certain symme- 
try for the average interface profile. This can be further 
elaborated by considering the continuum equation, and 
writing the interface field on one of the sides of the tri- 
angle as H (x, t) — m(t)x + SH(x, t). Inserting this into 
the QKPZ equation we obtain for SH 



asH(x, t) 
Ft 



v\7 z 8H 



\m{t)V5H + ^{V5Hf 



+v(x,H), (7) 



where the mean-field solution for H has been subtracted 
from the equation and the result is valid for intermediate 
times St <C t. We discover that the externally imposed 
slope (by the combination of F and H = at the bound- 
aries) changes the effective interface equation by creating 
a term linear in V5H, that will dominate over the KPZ- 
nonlinearity. The situation is thus exactly analogous to 
dcpinning transitions in the presence of an anisotropy, 
discussed by Kardar, Dhar, and Tang |jl5). This explains 
qualitatively why the avalanche behavior at criticality, 



at the SOC depinning transition, is different from that of 
the bulk dcpinning transition for the original equation. 
Regardless of whether the SOC ensemble differs from the 
translationally invariant one for ordinary SOC sandpiles, 
it is thus clear that in the case of SOC obtained with 
a similar recipe in QKPZ-systems the two ensembles do 
not share the same exponents. Should it be possible to 
write an activity-based description for the QKPZ univer- 
sality class, this implies that the field theory should have 
a different structure in the two cases, in contrast to the 
connection between the QEW and absorbing state phase 
transitions with a conserved field Jlq] , 

Above, we have concentrated on a cellular automaton 
for the QKPZ equation since as discussed it is related to 
a microscopic 'activity' picture, or a "QKPZ sandpile" . 
The slow-drive SOC-depinning limit can be taken also in 
the continuum version. Our investigations indicate that 
the same phenomenology of symmetry breaking persists 
in that case as well, though the actual exponents in the 
second phase before the critical slope is reached differ 
as Ah can now have arbitrary values and the avalanche 
properties change. The analysis of the continuum model 
is hampered by well-known numerical instabilities that 
develop once the smoothing effect of the Laplacian is re- 
duced - which is inevitable here due to the increasing av- 
erage slopes. We have also studied the two-dimensional 
case. While there, with periodic boundary conditions in 
one direction, and SOC ones in the other, are interesting 
details concerning the geometry of avalanche spreading, 
the basic physics is the same as in ID. 

In conclusion, in this paper we have analyzed the be- 
havior of interfaces in the QKPZ universality class in the 
presence of conditions that produce self-organized criti- 
cality. Spatial and temporal invariance are broken by the 
slow drive and the dissipativc pinning boundary condi- 
tions, giving rise to a macroscopic interface shape. The 
shape, in combination with the nonlinearity, results in 
the creation of a term linear in Sh in the effective equa- 
tion that governs the physics of avalanches. In contrast to 
usual SOC sandpiles we find that there are two avalanche 
regimes depending on the microscopic dynamics, where 
the second one fails to show critical scaling seen in the 
first. The ensemble in which the criticality appears is fun- 
damentally different from the usual dcpinning one and is 
related to the anisotropic QKPZ class Jl~5| ]. 

Partial support of OTKA T029985 is acknowledged. 
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